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SUMMARY 

A problem-adaptive solution procedure for Improving the reliability of finite element 
solutions to geometrically nonlinear shell-type problems Is presented. The strategy 
Incorporates automatic error detection and control and includes an Iterative procedure 
which utilizes the solution at one load step from one finite element model to obtain 
an equivalent solution at the same load step on a more refined model. Representative 
nonlinear shell problems are solved. 

INTRODUCTION 

Much of the research In adaptive finite element structural analysis has centered on the development 
of techniques for use In linear analysis. However, the need for an adaptive strategy Is even more 
Important In nonlinear analysis where a given finite element model may perform adequately for a 
certain range of loading (e.g., for the first several load steps) and become grossly Inadequate for 
another range of loading (e.g.. for the last several load steps). Recent work at the NASA Langley 
Research Center has focused on the development of an adaptive nonlinear analysis procedure for 
shell structures (McCleary (1]). This adaptive analysis procedure Integrates three primary compo- 
nents Into the nonlinear solution procedure: an automatic error detection strategy, an automatic 
error control strategy, and a reference state definition technique. Each of these components are 
described herein, and the use of the procedure Is demonstrated on two geometrically nonlinear shell 

problems. 

THE ADAPTIVE ANALYSIS PROCEDURE 

The primary components of the adaptive nonlinear solution procedure are Identified In the algorithm 
flowchart of Figure 1(a). The use of this algorithm to predict a nonlinear response curve Is 
Illustrated In Figure 1(b). The analyst defines an error tolerance which sets boundaries on the 
deviation of the calculated response from the converged response. At some load step, the Initial 
finite element discretization may predict a response which violates the prescribed error tolerance 
(e.g., see step 3 In Figure 1). At this point, a new finite element model Is generated and the 
solution procedure backs up to the previous load step (see steps 2 and 2' In the figure). Once an 
equivalent solution has been obtained at the previous load step for the new finite element mesh, 
the adaptive procedure proceeds until the response again violates the prescribed error tolerance 

(e.g., step 8' on the figure). 

Error Detection 

Many adaptive mesh refinement strategies employ a two tiered a posteriori error detection tec - 
nique. In this technique, a global error estimate Is computed as the summation of the local error 
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(a) (b) 

Figure l. Adaptive nonlinear solution procedure 


Indicators, and Is examined to determine If refinement of the finite element mesh Is necessary. 

When required, local error indicators are used to Identify the elements which need to be modified. 
Generally, the global error estimate is an accurate measure of the error In the finite element solu- 
tion. An alternate technique employs an error Indicator on an element-by-element basis to provide 
qualitative rather than quantitative Information about the solution, in this approach, local error 
Indicators are used only to guide mesh refinement. An estimate of the exact global error Is not 
calculated. 

This second approach has been adopted in the present study and may be viewed as element- 
independent In that element-specific information (e.g., element shape functions) is not accessed. 

The error Indicators are based on the use of stress resultants rather than the generalized nodal dis- 
placements. A smoothed stress resultant field, obtained by nodal averaging, is used to approximate 
the M exact M solution. As a consequence of using stress resultants, the error detection procedure 
Is Indirectly uncoupled from the large defiection/Iarge rotation formulation used to account for 
geometric nonlinearity. Specific knowledge of the form of the strain displacement relations (e.g., 
linear or nonlinear) or the nonlinear formulation (e.g., total Lagrange, corotatlonal) Is not required. 

As such, the procedure Is applicable to the geometrically nonlinear analysis of plate and shell 
structures. 

Two error Indicators have been employed In this study. The first is an adaptation of the energy norm 
of the error as proposed by Zienkiewlcz and Zhu [2J. The second indicator uses the discontinuities 
In the finite element stress resultant field In the form of a standard deviation In stress resultants. 

The error detection strategy developed by Zienkiewlcz and Zhu [2] for membrane problems and 
later extended to plate bending problems by Dow and Byrd [3], is based on the energy norm of the 
error. These element error indicators are defined herein by the energy norm of the error In stress 
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resultants within each element. Namely 

IMi = / (si- s;) T C(s; - s;)*). = / asJ <£ 1 AS k dn. (D 

where, for the k th element, SJ Is a vector of smoothed stress resultants (continuous across Interele- 
ment boundaries), S' k Is a vector of discrete stress resultants recovered from the finite element 
solution, and £„ Is the constitutive matrix that relates stress resultants to reference surface strain 
components. 

The energy norm of the error as defined In equation (1) Is an Integrated quantity and therefore 
requires either knowledge of the element shape functions or an element-independent generalization. 
The distribution of the difference between the smoothed and discrete values of the stress resultants, 
AS», over each element Is assumed to be described by Lagrange shape functions. The application 
of the energy norm of the error as an element error Indicator has at least two shortcomings. First, 
the evaluation of the Integral In equation (l) may become computationally expensive for large 
problems. Second, the error is In terms of the error In energy and thus provides only Indirect 
Information about the quality of the secondary solution (/.e., stress resultants). 

A new Indicator has been developed In an attempt to minimize the computational effort and to 
provide more direct Information about the quality of the stress resultant field. From statistics, the 
coefficient of variation of the stress resultants, V SR , Is defined as an error Indicator. In this work, 
V SR is defined at each node j In terms of the standard deviation In the stress resultants and the 
maximum absolute value of the smooth stress resultants (denoted as |S;| m „) as 

wSR _ 100 l T.?=i( S o, ~ S -i^ (2) 

_ is* | — , V 

The symbol n ( Is the number of elements connected to node j, and S„ Is a single stress resultant 
component. It should be noted that V® R Is defined for each stress resultant component. The 
elemental component value, vf R . Is taken to be the maximum nodal component value for a given 
element. The maximum coefficient of variation is used in order to ensure a conservative estimate 

CD 

of the error within an element as measured by V„ k . 

Error Control 

The development of an error control strategy requires the definition of refinement Indicators and 
the development of an automatic mesh refinement strategy. Refinement Indicators are used to 
Identify regions of the finite element mesh requiring an adjustment to their discretization and are 
formulated based on the previously defined error Indicators. The mesh refinement strategy, based 
on A-extenslon, determines precisely how these Identified regions are refined. 

The refinement Indicator based on the energy norm of the error is defined as In reference [2]. 
Elements needing refinement are Identified by using the local refinement Indicator ( k defined by 



The quantity ev» Is the projected elemental error assuming the current mesh Is an optimally refined 
mesh (/.e., the error uniformly distributed over the elements). The user-prescribed error tolerance 
is defined as a percentage of the global strain energy. 
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The refinement Indicator based on the coefficient of variation of stress resultants Is formulated 
by requiring that the analyst specify a maximum acceptable standard deviation in stress resultants 
as a fraction of the absolute maximum smoothed stress resultant. A single element refinement 
Indicator Is formed using the Euclidean norm, l.e., 



(4) 


where the error tolerance fi Is the fraction of the global maximum smoothed stress resultant. The 
use of an Euclidean norm In the definition of & ensures a conservative estimate of the refinement 
requirements for a given element. 


When the value of & exceeds unity In a given element, the element Is marked for refinement. A 
value of (h less than unity for a patch of elements Indicates that the patch should be less refined 
or undergo “fusion" (/.#., combining the elements In a patch of elements Into a patch containing 
fewer elements). Zlenklewlcz and Zhu (2] use an element sizing parameter, h/hk (equal to 1/& 
for 4-node elements and 1/V?T for 9-node elements), to predict multiple levels of refinement. The 
current Implementation contains no provisions for either multi-level refinement or the fusion of 
elements. 


An automatic mesh refinement strategy may take several forms and need not necessarily be fully 
adaptive. A fully adaptive strategy Is able to refine selectively. Automatic mesh refinement en- 
compasses any mesh refinement technique which can be performed automatically by the structural 
analysis software. However, automatic mesh refinement which Is not fully adaptive will generally 
tend to over-refine the finite element model. 

The automatic mesh refinement strategy employed In this study Is denoted as quasi-uniform mesh 
refinement since uniform refinement may be carried out within Individual, predefined regions of the 
finite element model. Quasi-uniform refinement operates best on finite element models which are 
easily parameterized so that It may be performed automatically by changing only the few parameters 
which define the finite element model. Within the framework of a general-purpose finite element 
code, quasi-uniform mesh refinement Is perhaps the most expedient approach in that It may be 
performed at a high level (e.g., using a command language), it does not require special data 
structures or software, and perhaps more Importantly, It always uses the precise definition of the 
structural geometry In the remeshlng process. 

Reference State Definition 

After a new finite element mesh has been defined, the nonlinear solution from the previous mesh 
needs to be defined for the new mesh. This solution on the new mesh corresponds to the reference 
state for restarting the nonlinear solution algorithm. Several factors must be considered In deciding 
the manner in which the reference state is defined. The reference state should be some Intermediate 
solution and not the Initial state (/.e., the solution procedure should not have to start from the 
first load step of the first finite element model). Computing the new reference state should not be 
computationally expensive, should not further complicate data structures, and should not require 
access to element-specific Information (such as shape functions). With these requirements In mind, 
a new approach has been developed In which Interpolation of displacement data Is not employed. 
The reference state is Instead defined using an Iterative solution recovery procedure [1]. 
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form applied displacement vector 
j — 1 

do for j < max_dJters 


uj 


form tangent stlfTness matrix based on uj_, 
solve uJ.j = f<_i using PCG method 
calculate error: • 

« = || u i_ t - ii / ii ui_ t ii 
If e < tot exit 

j = 3 + 1 
enddo 


= «fci 


Figure 2. Iterative solution recovery. 

The strategy for determining a consistent solution for a refined mesh Is summarized In Figure 2. 
This strategy uses the nonlinear solution at the previous load step (step « - 1) from the previous 
finite element model, to generate an applied displacement vector for the solution at the same load 
step on the new finite element model. Only those nodes which are new to the finite element model 
are left unconstrained during the Iteration cycle. An Initial estimate of the solution vector is formed 
by Initializing the displacements at the new nodes to zero. While these Initial displacements may 
be assigned values other than zero, doing so has not proven to alter significantly the convergence 
characteristics of the Iterative solution recovery procedure. Once formed, the Initial estimate is used 
to generate a tangent stlfTness matrix and serves as an Initial guess for a preconditioned conjugate 
gradient (PCG) Iterative solver. The solution Is compared to the solution at Iteration j - 1 
(l.e., end the solution recovery error Is measured by 
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(5) 


where Au<_i = u ti {f the error ' «£ ls ,ess than the spec,fied error to,erance ' a. converged 
solution at step i- 1 has been achieved. If the error Is greater than the specified tolerance, a new 
tangent stiffness matrix Is formed based on the j tk displacement solution. The procedure continues 
until a converged solution for load step i-lon the new finite element mesh is obtained. The 
algorithm then proceeds to load step i and obtains an Improved solution at load step i on the new 

finite element mesh. 

This technique Is particularly effective, and only efficient, when an iterative equation solver (e.g.. 
PCG technique) is used. In this case, the solution «C‘ Is used as an Initial solution for the Iterative 
solver at Iteration j (/.e.. as an Initial guess for u^). thereby Increasing the rate of convergence for 
the solver. While direct solvers may be used, they require a complete factorization of the tangent 
stlfTness matrix at each Iteration and will take roughly the same amount of time for each solution 
Iteration. An Iterative solver will generally require an ever decreasing number of Internal Iterations 
as the overall solution recovery process converges. 


NUMERICAL RESULTS 

The adaptive procedure described herein has been developed using the framework of the CSM 
Testbed [4j. Control of the analysis procedure Is implemented using the command language 
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Figure 3. Pear-shaped cylinder; geometry, properties, and response. 

capability of the CSM Testbed. Moreover, a stand-alone Fortran processor has been developed for 
evaluating the various error and refinement Indicators. A total Lagrangian formulation Is used for 
the nonlinear problem description, and the solution strategy Incorporates the Riks/Crlsfleld arc- 
length approach within a modified Newton-Raphson procedure to solve the system of nonlinear 
algebraic equations (see [1] for details). The refinement indicator, C*» displayed in the following 
sections Is analogous to the element sizing parameter previously discussed and Is defined as (t = 1/6 
for 4-node elements and as 6 = l/vt* for 9-node elements. 

In subsequent sections, the results of two nonlinear shell analyses are presented. Both problems 
involve structures which exhibit complex nonlinear behavior. The first problem, the collapse of a 
P«ar-shaped cylinder, demonstrates the use of the adaptive procedure on an isotropic, geometrically 
regular (rectangular elements are used throughout the domain) problem. Correlation with previously 
published results is shown. The second problem, the buckling of a cylindrical panel with a hole, 
demonstrates the use of the adaptive procedure on a composite, geometrically Irregular (skewed 
elements are required to model the central circular hole) problem. Correlations with test data and 
previously published analytical results are shown. 

Pear-Shaped Cylinder 

The pear-shaped cylinder shown In Figure 3 has been adopted by many researchers and the be- 
havior of this cylinder subject to a uniform end-shortening investigated (see Hartung and Bail 
[5]). The cylinder response becomes nonlinear at low values of applied end-shortening, and the 
normal deflections of the flat portions of the cylinder increase rapidly with increases in loading. 
The symmetry exhibited by the structure allows the analysis to be performed on one fourth of the 
cylinder. 

Quasi-uniform mesh refinement has been used for error control. The finite element mesh has been 
divided Into four groups of elements. Element group 1 contains the elements In the 45° arc at 
the top of the cylinder, element group 2 contains the elements in the flat segment between the 
two curved segments, element group 3 contains the elements in the 135° arc, and element group 
4 contains the elements In the bottom flat segment of the cylinder. Circumferential refinement 
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Figure 4. Refinement Indicators, C*. for Unite element models 1 through 5. 

Models A-E correspond to points A-E on Figure 3. 

within any Individual group of elements was unrestricted. Whenever all four element groups required 
refinement simultaneously, refinement through the depth of the cylinder was carried out in a uniform 
fashion (/.e. t for all elements in all groups). Command language procedures were used to direct 
the mesh refinement based on the original cylinder geometry. 

The nonlinear response of a series of five finite element models Is also shown In Figure 3. As 
shown In the figure, once the solution starts to deviate from the converged solution (603 nodes, 
132 9- node Assumed Natural— coordinate Strain (9ANS) elements [6J), mesh refinement occurs. 
The iterative solution recovery algorithm then defines the reference state for the restart of the 
nonlinear solution procedure. The collapse load is 2471 lbs., within 5% of the collapse load reported 
by Hartung and Ball [5]. 

For this series of models, refinement Is based on the energy norm of the error with a 10% error 
tolerance; each model Is composed of 9ANS elements. The refinement Indicators C* for the set of 
finite element models are shown In Figure 4. Initially, refinement was required only within the flat 
portions of the finite element model where significant nonlinear behavior was present. At higher 
values of load, the curved segments of the cylinder also needed refinement since these segments 
began to exhibit large normal deflections. Similar results are obtained using the coefficient of 
variation of stress resultants with an error tolerance of 10%. 

Composite Cylindrical Panel 

The postbucklfng response of axially compressed composite cylindrical panels with holes has been 
a subject of research for several years (e.g., Knight and Starnes [7]). This problem is characterized 
by large local deformations in the neighborhood of the hole which cause ply delaminations to occur 
In the postbuckling range. The panel analyzed herein is loaded with a uniform end-shortening 
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T30O6208 graphtte-epoxy taps 
16-ply quasMsotropic laminate 
15-inch radius 
14-inch squaro platform 
2-inch circular cutout 


Figure 5. Composite cylindrical panel: geometry, properties, and response 


and Is referred to as panel CP8 In reference [7j. The panel geometry and material and section 
properties are shown In Figure 5. 

To obtain a correct solution, the curvature of this panel must be preserved throughout the mesh 
refinement process. The exact panel geometry Is difficult to represent because of the presence of a 
cylindrical surface with a curved Interior boundary. In general, the finite element mesh will require 
many distorted elements In order to adequately represent the panel geometry. A uniform mesh 
refinement strategy was used to accomplish mesh refinement. Once any single element was flagged 
for refinement, uniform mesh refinement was carried out by alternately adding spokes of nodes and 
rings of elements. A command language procedure was used to direct this mesh refinement. 

The analyses presented in this section were halted soon after buckling occurred due to the in- 
adequacy of the material model after buckling (as reflected in the experimental results). The 
converged solution shown In Figure 5 (which plots the normalized load versus end-shortening re- 
sponse) was obtained using a finite element mesh with 1600 nodes and 384 9ANS elements. While 
the predicted buckling load is below the experimental data, it Is consistent with the predicted 
buckling load given In reference (7]. in addition, the effects of geometric Imperfection were not 
included In the analysis which may account for some of the differences between the experimental 
and analytical results. 

A series of four finite element meshes was generated using a 15% error tolerance on the coefficient 
of variation of stress resultants. The model used 4-node Assumed Natural-coordinate Strain 
(4ANS) elements [6]. The resulting response curve is also shown In Figure 5. Again the solution 
recovery algorithm had no difficulty in transitioning from one Unite element mesh to another. Three 
finite element meshes were automatically generated before the error Indicators computed at the 
first load step indicated that an acceptable solution had been obtained. The refinement Indicators 
for this series of finite element meshes are shown In Figure 6 in a planform view. As the model Is 
refined, the elements requiring refinement move along the diagonals of the panel toward the edge 



Refinement Indicators, <*, for finite element models 1 through 4. 
Models A-D correspond to points A-D on Figure 5. 


of the hole where large gradients occur. 

Computati onal Effort 

One of the goals of the work described thus far Is to show that the adaptive nonlinear solution 
strategy presented herein can be efficient. Satisfying this goal requires that the computational 
time required for performing a nonlinear analysis using the adaptive procedure as outlined should 
not substantially exceed the computational time required for a nonlinear analysis performed using 
a single refined finite element model. Timing data for both the pear-shaped cylinder and the 
composite cylindrical panel are summarized In Table 1. For each case, the single refined finite 
element model was selected to be the final finite element model generated by the adaptive solution 
procedure. This single refined model was defined only after the adaptive analysis was performed. 
In general, multiple models and analyses would have been required to obtain a reliable solution. 

in the case of the pear-shaped cylinder, the analysis performed using the adaptive procedure 
required only approximately 9% more CPU time than the analysis performed using the single refined 
finite element model. In the case of the composite cylindrical panel, the analysis performed using 
the adaptive procedure required approximately 19% more CPU time than the analysis performed 
using a single rellned model. This difference In the additional percentage of CPU time required Is 
attributable to the fact that all of the mesh refinement occurs within the first two load steps of 
the analysis of the composite panel. While the analysis performed using the single refined model 
required less CPU time for both problems, the reliability of the solutions can only be assessed 
by a complete re— solution on a different mesh. In both cases, the total CPU time required for 
computing two complete nonlinear solutions would be greater than the CPU time required to 
perform the single adaptive analysis. 

The results shown In the table are very encouraging In that the adaptive analysis Is at least 
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Tabic 1. Timing Date for Adaptive and Single Non linear Analyse* 



Pear-Shaped Cylinder 

Composite Cylindrical Panel 

Function 

Adaptive 

Analysis 

Timet 

Single 

Refined Model 
Timet 

Adaptive 

Analysis 

Timet 

Single 

Reflned Model 
Timet 

Model Generation 
Form and Factor 

100.6 

4.2 

81.7 

34.5 

Stiffness 

2561.2 

2497.5 

3254.9 

3062.0 

Matrix Operations 

967.0 

949.7 

584.5 

575.0 

Error Calculation 

117.0 

- 

435.3 

- 

Total 

3746.1 

3452.8 

4358.6 

3672.5 


t In CPU seconds. Calculations performed on a Convex C220 minisupercomputer. 


competitive with the analysis performed on the single refined finite element model, if a fully 
adaptive (/.e. f selective) mesh refinement strategy was employed, at least two Improvements in the 
efficiency of the procedure would be possible. First, the model generation phase of the analysis 
would be more efficient. The current Implementation requires a complete regeneration of each 
model rather than an Incremental change In the previous model. Second, improvements could 
be made In the Iterative solution recovery procedure by partitioning the tangent stlfTness matrix 
and only reforming the elemental matrices for the "new" elements. The current Implementation 
requires the complete regeneration of the tangent stiffness matrix for each solution Iteration. 

CONCLUDING REMARKS 

A problem-adaptive nonlinear solution procedure has been described. This procedure Incorporates 
automatic error detection and control Into a modified Newton-Raphson nonlinear solution strategy. 
A technique for defining a reference state, called Iterative solution recovery, has been developed. 
This technique has been shown to be an effective means of making the transition from one finite 
element discretization to another at a given load level. The use of two error indicators, the energy 
norm of the error and the coefficient of variation of stress resultants, as effective guides for mesh 
refinement, has been demonstrated. Two geometrically nonlinear shell problems have been solved 
and good correlation with published results has been presented. 
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